superplot <- function (results, start, end, rows, columns, indiv = 1, lmer = 0, var.slope = 0, fiftydots = 0, rich = 7, med = 35, poor = 24, fiftylines = 0, survey="", predictor="Individual Income", lowbound=-2,hibound=2, inc.bins=1, xlabels="all", ylabels="all", sf=4.5) 
{
    #i.result=round(br.state[[13]]$summary,2)

    #windows(height=8,width=6.5)
    par (mfrow = c(rows,columns))
    par(mar=c(3.75, 4.25, 3, 1)) # bottom, left, top, right
    
    for (j in start:end) {
        if (ylabels[1] == "all") {        
            if (indiv == 0) { ylabel = "Rep Vote Share" } else { ylabel = "Probability Voting Rep" }
        } else {
            if (j %in% ylabels) { ylabel = "Probability Voting Rep" } else { ylabel = "" }        
        }
        
        if (xlabels[1]=="all") { 
            xl = predictor 
        } else {
            if (j %in% xlabels) { xl = predictor } else { xl = "" }
        }
        
        plot(c(lowbound,hibound), c(.1,.9), xaxt="n", yaxt="n", xlab=xl, ylab=ylabel, cex.lab=1.5, mgp=c(2.5,.5,0), type="n")
        
        axis(side=1,seq(lowbound,hibound,(hibound-lowbound)/4),cex.axis=1.5)
        axis(side=2,seq(.25,.75,.25),cex.axis=1.5)
        #axis(side=2,seq(.12,.9,.2),tick=T,labels=F)
        title(paste(survey,main=years[j],sep=""), cex.main=1.75)
        #range(st.repshare[st.year==years[j]])
        
        superp = matrix(nrow=50,ncol=9)
        rownames(superp)=state.name
        colnames(superp)=c("Intercept","Slope","W Mean","Pred","1","2","3","4","5")     
        
        #attach(cps.data[[j]])
        cps.stats=cps.stats.all[[j]]
        
        if (lmer == 1 ) {
           
            constant = 0; slope = 0; # common constant
            #var.int.old = results[[j]][,1]
            #var.slope.old = results[[j]][,2]
                        
            #if (j == 5) {
                statelist = as.integer(rownames(results[[j]]))
                #statelist = seq(1:50)
                #st.income = log(st.inc10k[st.year==years[j]])[statelist]
                #st.income = st.income[as.integer(statelist)]
                #var.int <- results[[j]][,1] + results[[j]][,3]*st.income
                #var.slope <- results[[j]][,2] + results[[j]][,4]*st.income
                var.int <- results[[j]][,1]
                var.slope <- results[[j]][,2]

                #names(var.int)[8:length(var.int)]=as.integer(names(var.int))[8:length(var.int)]-1
                #names(var.slope)[8:length(var.slope)]=as.integer(names(var.slope))[8:length(var.slope)]-1
                #statelist = names(var.int)
                           
            #} else {
                #st.income = cps.stats[,4][1:49] # No Wyoming
                #statelist = as.integer(rownames(results[[j]]))
                #statelist = seq(1:50)
                #st.income = log(st.inc10k[st.year==years[j]])[statelist]
                #var.int <- results[[j]][,1] + results[[j]][,3]*st.income
                #var.slope <- results[[j]][,2] + results[[j]][,4]*st.income
                #var.int <- results[[j]][,1]
                #var.slope <- results[[j]][,2]
                #names(var.int)=seq(1:50); names(var.slope)=seq(1:50)
                #statelist = names(var.int)
 
            #}

                        
            for (i in 1:50) {
                if (sum(statelist==i) > 0) {
                    superp[i,"Intercept"] = var.int[statelist==i]
                    superp[i,"Slope"] = var.slope[statelist==i]
                    superp[i,"W Mean"] = cps.stats[i,"Weighted Count"]
                    superp[i,"Pred"] = round(invlogit(superp[i,"Slope"]*superp[i,"W Mean"]+superp[i,"Intercept"]),3)
                    superp[i,5:9] = cps.stats[i,5:9]
                }
            }
            #superp[9:49,1:2]=superp[10:50,1:2]
 
            # end of lmer
            } else {
            if (indiv == 1) {
                if (is.vector(results[[j]])) { means = results[[j]] } else { means = results[[j]][,"50%"] }
                constant = 0; #constant=means["beta.0.adj"]
                slope = 0; #slope=means["b.inc.0.adj"]
                offset = 50
            } else {
                means = results[[j]][,"mean"]
                constant=means["constant.adj"]
                if (var.slope == 0) { slope = results[[j]][52] }
                offset = 0
            }

            for (i in 1:50) {
                if (indiv == 0) { superp[i,"Intercept"] = constant } else { superp[i,"Intercept"] = means[i]+constant }
                if (var.slope == 0) { superp[i,2] = mean(results[[9]][,"50%"][51:100]) } else { superp[i,2] = means[offset+i]+slope}
                superp[i,3] = cps.stats[i,"Weighted Count"]
                if (indiv == 0) {
                    superp[i,"Pred"] = round(slope*superp[i,"W Mean"]+superp[i,"Intercept"],3)
                } else {
                    superp[i,"Pred"] = round(invlogit(superp[i,"Slope"]*superp[i,"W Mean"]+superp[i,"Intercept"]),3)
                }
                #superp[i,5:9] = round(count[fips.cps==s2f[i]]/cps.stats[i,"Sum Count"],2)
                superp[i,5:9] = cps.stats[i,5:9]
                }
        # end of non-lmer
        }
 
        if (fiftydots==1) {
            for (i in 1:50) {
                points(superp[i,3],superp[i,"Pred"],cex=1.25,pch=19,col="gray")
            }
            fit <- lm (superp[,"Pred"] ~ superp[,"W Mean"])
            abline (fit$coeff[1],fit$coeff[2],col="gray")
        }
        
        abline (.5,0,col="gray")


        if (fiftylines == 0) {
            if (indiv == 0) {
                abline (superp[rich,1], superp[rich,2], col="blue") # rich
                abline (superp[med,1], superp[med,2], col="green") # medium
                abline (superp[poor,1], superp[poor,2], col="red") # poor
            } else {
                curve (invlogit (superp[rich,1] + superp[rich,2]*x), lowbound, hibound, add=T) #, col="blue") # lty = 1
                curve (invlogit (superp[med,1] + superp[med,2]*x), lowbound, hibound, add=T) #, col="green") # lty = 2
                curve (invlogit (superp[poor,1] + superp[poor,2]*x), lowbound, hibound, add=T) #, col="red") # lty = 3
            }
        } else {
        print ("fiftylines!\n")
        for (i in 1:50) {
            curve (invlogit (superp[i,1] + superp[i,2]*x), lowbound, hibound, add=T, lwd=0.5, col="blue")
            }    
        }
        
                                
        if (inc.bins == 0) {
            #for (i in c(rich,med,poor)) {
                #if (i==rich) { offset = .03 } else { offset = -.01 }
                #text(1.45,invlogit(superp[i,2]*2+superp[i,1])+offset, state.name[i],cex=1.35)
            inc.cat = c(-1,-.6,-.2,.2,.6)        
            st.label = .4
        } else {
            inc.cat = c(-2,-1,0,1,2)
            st.label = 1.35
        }
        
            #sf = 4.5
    
            # draw circles for income quintiles
            col.inc = "black"
            if (indiv == 0) {
                for (i in c(rich,med,poor)) {
                    points(superp[i,"W Mean"],superp[i,"Pred"],cex=sf,pch=19)
                    points(inc.cat[1],slope*inc.cat[1]+superp[i,1],cex=sqrt(superp[i,"1"])*sf,col=col.inc,pch=21)
                    points(inc.cat[2],slope*inc.cat[2]+superp[i,1],cex=sqrt(superp[i,"2"])*sf,col=col.inc,pch=21)
                    points(inc.cat[3],slope*inc.cat[3]+superp[i,1],cex=sqrt(superp[i,"3"])*sf,col=col.inc,pch=21)
                    points(inc.cat[4],slope*inc.cat[4]+superp[i,1],cex=sqrt(superp[i,"4"])*sf,col=col.inc,pch=21)
                    points(inc.cat[5],slope*inc.cat[5]+superp[i,1],cex=sqrt(superp[i,"5"])*sf,col=col.inc,pch=21)
                    }
            } else {
                for (i in c(rich,med,poor)) {
                    points(superp[i,"W Mean"],superp[i,"Pred"],cex=sf,pch=19)
                    points(inc.cat[1],invlogit(superp[i,2]*inc.cat[1]+superp[i,1]),cex=sqrt(superp[i,"1"])*sf,col=col.inc,pch=21)
                    points(inc.cat[2],invlogit(superp[i,2]*inc.cat[2]+superp[i,1]),cex=sqrt(superp[i,"2"])*sf,col=col.inc,pch=21)
                    points(inc.cat[3],invlogit(superp[i,2]*inc.cat[3]+superp[i,1]),cex=sqrt(superp[i,"3"])*sf,col=col.inc,pch=21)
                    points(inc.cat[4],invlogit(superp[i,2]*inc.cat[4]+superp[i,1]),cex=sqrt(superp[i,"4"])*sf,col=col.inc,pch=21)
                    points(inc.cat[5],invlogit(superp[i,2]*inc.cat[5]+superp[i,1]),cex=sqrt(superp[i,"5"])*sf,col=col.inc,pch=21)
                    #if (i==rich) { offset = .1 } else { offset = -.01 }
                    text(st.label,invlogit(superp[i,2]*2+superp[i,1] + .1), state.name[i],cex=1.25)
                    }
            }
                    
            #legend (1,.7,c(state.abb[rich],state.abb[med],state.abb[poor]),col=c("green","brown","orange"),lty=1, lwd=3, bg="gray85")
            
            lines(lowess(c(superp[rich,"W Mean"],superp[med,"W Mean"],superp[poor,"W Mean"]), c(superp[rich,"Pred"],superp[med,"Pred"],superp[poor,"Pred"]) ))
            
        }

    return (superp)
}
